/* lambert.cpp */

#include <dodo_plot/map/lambert.h>
#include <dodo_plot/kernel/position.h>
#include <algorithm>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

namespace dodo_plot
{
using namespace std;

Point Lambert::operator() (double lon, double lat) const
{
	double rho_ = f_ /
		pow (
			tan(M_PI / 4 + lat * M_PI / 180.0 / 2),
			n_
		    );

	double x = rho_ * sin( n_ * (lon * M_PI / 180.0 - ref_lon_) );
	// Note: this y increase as latitude increase.
	// It's opposite against the Cairo axis.
	double y = rho_0_ - rho_ * cos( n_ * (lon * M_PI / 180.0 - ref_lon_));
	
	return pos_proj_(x, y);
}

Point Lambert::operator() (const Point& lonlat) const
{
	return operator() (lonlat.x, lonlat.y);
}

Point Lambert::map_proj_only(double lon, double lat) const
{
	// This is only for initially the Pos_Projection
	double rho_ = f_ /
		pow (
			tan(M_PI / 4 + lat * M_PI / 180.0 / 2),
			n_
		    );

	double x = rho_ * sin( n_ * (lon * M_PI / 180.0 - ref_lon_) );
	// Note: this y increase as latitude increase.
	// It's opposite against the Cairo axis.
	double y = rho_0_ - rho_ * cos( n_ * (lon * M_PI / 180.0 - ref_lon_));

	return Point(x, y);
}

Lambert::Lambert(double standard_parallel_1,
                double standard_parallel_2,
                double ref_lon,
                double ref_lat,
		const Point& corner_1,
		const Point& corner_2
               ) :
	std_par_1_(standard_parallel_1 * M_PI / 180.0),
	std_par_2_(standard_parallel_2 * M_PI / 180.0),
	ref_lon_(ref_lon * M_PI / 180.0),
	ref_lat_(ref_lat * M_PI / 180.0)
{
	init();
	Point corner_xy_1 = map_proj_only(corner_1.x, corner_1.y);
	Point corner_xy_2 = map_proj_only(corner_2.x, corner_2.y);
	
	double ll = std::min(corner_xy_1.x, corner_xy_2.x);
	double rr = std::max(corner_xy_1.x, corner_xy_2.x);
	// notice the lambert projection y direction is opposite against Cairo.
	double tt = std::max(corner_xy_1.y, corner_xy_2.y);
	double bb = std::min(corner_xy_1.y, corner_xy_2.y);
	pos_proj_ = Pos_Projection(Position(0, 1.0, 0, 1.0), Position(ll, rr, tt, bb));
}

void Lambert::init()
{
	n_ = log( cos(std_par_1_) / cos(std_par_2_) ) /
		( 
			log (
				tan(M_PI / 4 + std_par_2_ / 2) /
				tan(M_PI / 4 + std_par_1_ / 2)
			)
		);

	f_ = cos(std_par_1_) /
		pow(
			tan(M_PI / 4 + std_par_1_ / 2),
			n_
		   ) /
		n_;

	rho_0_ = f_ /
		pow(
			tan(M_PI / 4 + ref_lat_ / 2),
			n_
		   );
}

}
